/*
 * Cholesky.cpp
 *
 *  Created on: 2014年6月27日
 *      Author: fangxia722
 */

#include "Cholesky.h"
#include "Matrix.h"
#include "Error.h"
#include <string.h>
#include"BandMatrix.h"
/**
 * 使用基于外积的 Cholesky分解
 * PAP(T)=GG（T）
 * 可以适应半正定矩阵（包括数值的半正定矩阵），若非正定矩阵抛出异常
 * 包括选主元计算
 */
Matrix* Cholesky::decom(Matrix& A) {
	if (!A.isSymmetric()) {
		return NULL;
	}
	Matrix* result = A.clone();
	if (type == 0) {
		for (int k = 0; k < result->getColL(); k++) {
			if (result->get(k, k) > 0) {
				result->set(k, k, sqrt(result->get(k, k)));
				for (int i = k + 1; i < result->getRowL(); i++) {
					result->set(i, k, result->get(i, k) / result->get(k, k));
				}
				for (int j = k + 1; j < result->getColL(); j++) {
					for (int i = j; i < result->getRowL(); i++) {
						result->set(i, j,
								result->get(i, j)
										- result->get(i, k)
												* result->get(j, k));
					}
				}
			}
		}
	} else {
		if (p != NULL) {
			delete p;
		}
		p = new int[A.getColL()];
		memset(p, 0, sizeof(int) * A.getColL());
		//选主元所用变量
		int r = 0;
		for (int k = 0; k < result->getColL(); k++) {
			//寻找合适的q(k≤q≤n)使得A(q,q)=max|A(k.k),⋯A(n,n)|
			int q = k;
			for (int qq = k + 1; qq < result->getColL(); qq++) {
				if (result->get(qq, qq) > result->get(q, q)) {
					q = qq;
				}
			}
			p[k] = q;
			if (result->get(q, q) > 0) {
				r = r + 1;
				if (k != q) {
					result->exchangeRow(k, q);
					result->exchangeCol(k, q);
				}
				//正常的分解过程
				result->set(k, k, sqrt(result->get(k, k)));
				for (int i = k + 1; i < result->getColL(); i++) {
					result->set(i, k, result->get(i, k) / result->get(k, k));
				}
				for (int j = k + 1; j < result->getColL(); j++) {
					for (int m = j; m < result->getColL(); m++) {
						result->set(m, j,
								result->get(m, j)
										- result->get(m, k)
												* result->get(j, k));
					}
				}
			} else if (result->get(q, q) < 0) {
				throw SYSTEM_NOT_SUPPORT;
			}
		}
	}
	return result;
}

